Mean-field vs. stochastic models for transcriptional regulation 
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We introduce a minimal model description for the dynamics of transcriptional regulatory net- 
works. It is studied within a mean-field approximation, i.e., by deterministic ode's representing the 
reaction kinetics, and by stochastic simulations employing the Gillespie algorithm. We elucidate 
the different results both approaches can deliver, depending on the network under study, and in 
particular depending on the level of detail retained in the respective description. Two examples are 
addressed in detail: the repressilator, a transcriptional clock based on a three-gene network realized 
experimentally in E. coli, and a bistable two-gene circuit under external driving, a transcriptional 
network motif recently proposed to play a role in cellular development. 

PACS numbers: 87.18.Cf, 87.10.Ed, 87.10.Mn 



I. INTRODUCTION 

Mathematical models for the dynamics of transcrip- 
tional regulation are traditionally formulated either in 
terms of ordinary differential equations 0,01, or by 
purely stochastic models, based on Master equations [3j 
or by using the Gillespie algorithm 0]. Both the deter- 
ministic and stochastic descriptions average out spatial 
degrees of freedom and hence are more similar to each 
other than is often acknowledged. In recent years, a dis- 
cussion has started on the effect of stochasticity on gene 
regulatory processes; exemplary studies are [U IS H @ • 
Indeed, already the fact that molecules involved in reg- 
ulatory processes often exist only in small copy numbers 
can be significant for the dynamics of a given regulatory 
circuit, and stochastic effects like bursting may have an 
important role for cellular function 0. 

Models of regulatory dynamics suffer also from another 
problem which is the lack of precise knowledge of reac- 
tion rates. Building dynamic models for a large num- 
ber of network elements can induce further arbitrariness 
due to a lack of detailed knowledge of the interaction 
mechanisms involved. Approaches that aim to describe 
larger networks are often deliberately reductionist to be- 
come computationally tractable (see, e.g., build- 
ing on p ioneering work by Glass, Kauffman and Thomas 
111 lQ)) an d the result of such computations can then 
only be called "qualitative". The effect of these reduc- 
tion schemes, which within a physics-based notion could 
also be subsumed under the notion of "coarse-graining" , 
therefore often lacks clarity as to what effect the approxi- 
mations/simplifications have, since a general systematics 
is not available (an exemplary discussion of this issue can 
be found in [Hj]). 

In this paper we address the question of what effect 
such a reduction scheme has on the dynamics of a given 
regulatory network in a systematic way. For this we start 
from a minimal model description for transcriptional reg- 
ulatory networks which coarse grains as many regulatory 
layers as possible (although they could of course be added 



back in later). We note that this modeling philosopy is 
in contrast to the usual way models of transcriptional 
regulation are built in which first all available biochem- 
ical detail is considered and then reduced by way of ap- 
proximation (as, e.g., in [Til [l5j and many other simi- 
lar examples). We then formulate both a deterministic 
(mean-field) version and a stochastic version of the tran- 
scriptional dynamics. This approach allows us to study 
the dynamics of basically all fundamental classes of tran- 
scriptional networks relevant for prokaryotic organisms, 
although wc only look at few-gene networks in detail here. 

The outline of the paper is as follows. We first develop 
the kinetic reactions involved in transcriptional regula- 
tion. Subsequently, we formulate the corresponding de- 
terministic and stochastic versions of the dynamics. A 
separate section of the paper is devoted to the applica- 
tion of both schemes to commonly encountered regula- 
tory motifs Two examples are presented in more 
detail since they display richer structure: the repressila- 
tor, a three-gene network of inhibiting gates which acts 
as a genetic clock, previously realized experimentally in 
E. coli 17], and a regulatory motif with multiple inputs 
which was recently proposed to be relevant for regulatory 
processes in development [llj]. For all these systems, we 
compare the results of the deterministic calculations and 
their stochastic counterparts and evaluate the role dif- 
ferent regulatory mechanisms play for the observed out- 



II. THE GENE GATE MODEL 
A. The transcriptional reactions 

Our minimal model for transcriptional regulation con- 
sists in the definition of a computational element for each 
regulatory element (i.e., transcribing gene), which we call 
a gene gate. The basic possible types of gene gates are 
sketched in Figure 1. Each gene gate is defined via its 
reaction kinetics. The 'null gate' in Figure 1.1 is a gene 
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Figure 1: (Color online) The four basic types of gene gates: 
1) The null gate (a gate without control input); 2) The neg 
gate (repression of transcription); 3) The pos gate: activation 
of expression; 4) the posneg gate: a multi-input gate with one 
activating and one repressing input. 



and translation are lumped together. 

Likewise we can model the activation of a gene upon 
binding of a transcription factor; Figure 1.3 shows the 
'pos gate'. The binding reaction is identical, but the 
gene in state G' now behaves according to 



G' G + B 



(5) 



where the rate 77 > e, i.e. the transcription/translation 
rate upon activation is larger than the basal rate. This 
is the pos (A; S)-gate. 

Finally, Figure 1.4 shows a gate with multiple regula- 
tions which is in fact a commonly encountered situation, 
see, e.g., the E. coli network of transcriptional interac- 
tions reconstructed in [l9[. For the posneg(A, C; -B)-gate 
we have to consider three gene states, G, G', and G" with 
the reactions 



G + G 



G' + C 



(0) 



in a state G which produces a protein output B at a rate 
e, hence the kinetic reaction is written as 



A + G 



G" 



G 



G + B. 



(1) 



The protein output can be degraded according to the 
reaction 



B 



0. 



(2) 



In an abbreviating notation we call this gate element 
null(0;B) where inputs and outputs are separated by 
the semicolon. 

In the next step we add a regulatory input to the null 
gate. Figure 1.2 shows the resulting 'neg gate' in which a 
transcription factor A inhibits the production of protein 
B upon binding. This is represented by the reaction 



A + G^ r G' + A. 



(3) 



This reaction corresponds to the formation of a tran- 
scription factor-DNA complex with zero lifetime; such 
an intermediate with a finite lifetime can of course be in- 
troduced as well but is not necessary for a minimal model 
of gene networks. 

After this interaction, the gene in state G' is blocked in 
transcription/translation. In order to allow transcription 
again the gate has to relax from its blocked state to its 
original transcribing state at a rate 77, 



G 1 



G 



(4) 



to the state G in which transcription at a basal rate e can 
occur. We call this gate the neg(A; S)-gate. The relax- 
ation process from G' to G models the fact that a gene 
generally is not transcribed immediately after the break- 
up of a transcription factor-DNA complex; also note that 
within our minimal model of the gene gate, transcription 



and the correponding relaxation reactions 
G —>rii G 



G" 



G + B. 



(7) 



(8) 



(9) 



It is clear from this scheme that for each additional regu- 
latory function, a binding transcription factor and a cor- 
responding gene state have to be introduced. 

Our minimal model obviously leaves out a number of 
regulatory levels such as 

• complexation of transcription factors; 

• formation of the DNA-transcription factor com- 
plex; 

• DNA transcription and RNA translation are 
lumped together. 

These regulatory mechanisms can, of course, be added to 
the list of reactions given above, and we will come back 
to this issue in the course of this paper. 



B. The mean-field equations 

Having listed the transcriptional reactions we now de- 
fine a continuum description based on ordinary differen- 
tial equations for the concentration of genes and proteins. 
We assume that the cell population can be considered as 
a 'soup' containing the proteins as well as N copies of 
the gene G. We denote normalized concentrations by 
small letters g = [G]/N, b= [B]/N with [G] = #G/V 
(likewise for [B]) and keep the previous symbols for the 
kinetic constants (i.e., we include dependencies on cell 
volume V and gene copy number N where necessary; the 
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difference to the kinetic reactions should be evident from 
the context). The two reactions of the null gate are then 
summarized by the ode 



eg - 76 . 



(10) 



For the regulated genes, an equation for g has to be 
added. Since the iV gene gates present in our cell model 
have to be either in state G or G' , one has the conser- 
vation law [G] + [G'\ = N. From the normalization we 
have g + g' = 1, and hence the neg-gate is described by 
the two odes, eq. tjTt)]) . and 



9 = Vf/ - rga = 77(1 - (1 + va)g) . 



(11) 



where the conservation condition has been used, and v — 
r/rj. 

The pos-gate (Figure 1.3) is governed by the ode's 
eq. flTTJ) and 



b = £.9 + V9' - J b = V - (V - e)s - l b ■ 



(12) 



Finally, we consider the case of multiple regulations of 
a single gene, the simplest multi-input gate, the posneg- 
gate of Figure 1.4 with the three gene states, G, G' and 
G", modifying the conservation condition to g + g' +g" — 
1. We can build up the gate reaction kinetics as before 
and obtain the system of ode's 



and 



b = eb + rj 2 g" - ■yb . 



-Vi9 +ngc, 



9 



-V29 +r 2 ga : 



(13) 



(14) 



(15) 



hence one has for g the equation g = —(g' + g") which 
follows from the conservation of gene states. 

At this point we stress that we have only considered 
the case of binding of a single protein A. In general, the 
binding of proteins is rather by multi-protein complexes 
(dimers or higher) , which is one way to give rise to a Hill 
coefficient h when the complexation reaction is consid- 
ered an equilibrium ( "fast" ) reaction [20] . We could take 
this into account in our model by adding a correspond- 
ing complexation reaction in the reaction scheme. To be 
practical we here directly modify the ode equation of the 
gene by replacing a by a h with h > 1 to cover this more 
general case; in what follows, we consider h as a conti- 
nously variable parameter. It is well-known that a Hill 
exponent > 1 is essential for the dynamic behaviour of 
simple gene circuits (2lj . 

For the stochastic simulations we employ the Gille- 
spie algorithm which is equivalent to the Chemical Mas- 
ter equation Q. We combine the Gillespie method with 
the stochastic 7r-calculus, a process algebra originating in 
theoretical computer science [H, [H, [H HI [27J . For 
a brief introduction into the main ideas of the calculus, 



see Appendix A. 



III. EXAMPLES 



Basic circuits 



We first discuss the elementary gene circuits that can 
be built from the above constructs. All simple transcrip- 
tional networks are either circular, linear or mixed cir- 
cuits, see Figure 2. The archetypal loops are the autoin- 
hibitory and autoactivatory loops. The autoinhibitory 
loop neg(a;a) is shown in Fig 2.1. The ode's governing 
its dynamics are 



a = eg-ya : 



and 



.9 = V9 



rga = ?/(l - (1 + va h )g) . 



(16) 



(17) 



The natural first task is to look at nullclines and fixed- 
points. The nullcline of g is determined by 



1 



9 = 



1 



(18) 



If we have g/rj w and v finite we can keep the circuit 
near the nullcline of g. Inserting the nullcline condition 
into the equation for a we find 



7a, 



(19) 



which is the common form of the Hill-type equation 
used in nonlinear dynamics descriptions of gene networks. 
This turns out to be a general feature of the gene gate 
approach: near the nullclines of the gene gate states, 
g w g' w .... w 0, the circuit dynamics reduces to that of 



£3- r^Q_ ITTrT 
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Figure 2: (Color online) The two main classes of simple cir- 
cuits: circular (1) and linear (2). Shown are only the repres- 
sive circuits; activatory circuits and mixtures of both types 
can be built in a similar fashion. Circuits shown in (1): the 
autoinhibitive circuit, a bistable switch, the repressilator. Cir- 
cuits in (2): a linear array and a linear array with a head 
feedback: hence a mixture of a circular and a linear circuit. 
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the standard Hill equations. This feature has an imme- 
diate consequence for the fixed points. The nuclline of a 
is given by 



1 



la, 



(20) 



where the result for g has been used, and we thus find 
the standard fixed-point condition of the Hill equation 
for a. Since the left-hand side is a hyperbolic function in 
a, and the right-hand side is a linear function there is a 
unique fixed-point of the circuit. 

The argument can be repeated for the autoactivatory 
loop pos(a; a) with the result 



n 



v 



l 



"fa 



1 



(21) 



which is the typical sigmoidal form of the activatory cir- 
cuit. Therefore, we again find that the fixed-points are 
given by a conditions akin to the standard Hill-type equa- 
tions, which for h > 1 gives rise to three fixed-points. 

The stability of the fixed-points in the gene networks 
is not affected by the presence of the genes. We illustrate 
this for the bistable circuit composed of two neg-gtaes, 
neg(a\ b)\neg(b\ a), where the symbol | denotes the com- 
position of two gates, see Figure 2.1. The equations of 
th circuit read as 



and 



a = eg a - ja 



g a = 77(1 - (1 + vb h )g a ) 



(22) 



(23) 



and likewise for a <-» b. As is well known [2l|, the non- 
linearity due to the Hill coefficient is needed for the sys- 
tem in order to display the fixed-point structure of the 
bistable switch; for a value of h — 1 as in our basic ver- 
sion of the gene gate model this is not the case. The 
stability of the fixed-points follows from the eigenvalues 
of the matrix 



r fp = 



-1 


e 











-x 


-£ 











-7 





-£ 








-x 



(24) 



Taking the root of this equation, one finds four real eigen- 
values, two of which are negative, and two positive. The 
picture that emerges therefore is the usual instability in 
the space of protein concentrations a\ , a2, while the genes 
do not contribute. 

We close this subsection by commenting on results 
from the stochastic simulations. The basic loop- and 
linear circuits (negative, positive) show fixed-point be- 
haviour similar to their deterministic counterparts (24J. 
For the bistable switch there is a notable difference: as 
was recently shown based on a Master equation approach 
the stochastic dynamics of the bistable switch without co- 
operativity (h = 1) displays both bistability and switch- 
ing [28j | . This behaviour is easily reproduced with our 
Gillespie approach, see Figure 3. 
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Figure 3: (Color online) Switching in the stochastic bistable 
circuit without cooperativity. Simulation parameters are: r 
=1, e — 0.4, 77= 0.2, 7 = 5 -10 -3 . The insert indicates out- 
put on the 7r-calculus channels a!, 6!, equivalent to protein 
numbers. 



Before moving on to richer examples, we draw a brief 
intermediate conclusion for the gene gate model: 



with 



X EE 77(1 + ^) , £ 



rha h 

1 + V( 



(25) 



Note that we are looking here at the stability of the sym- 
metric fixed-point for which \i = X2, £1 = £2- For the 
bistable switch, this is the unstable fixed-point interven- 
ing between the two stable fixed-points, and its eigenval- 
ues follow from the characteristic polynomial to T/p, 



( 7 + A) 2 (A + x) 2 -(£0 2 



(26) 



• if the deterministic gene circuit has a unique stable 
fixed-point, the genes are 'irrelevant' variables in 
the sense that they do not alter the location of the 
fixed point. They do, however, affect the transient 
dynamics (see below); 



• the deterministic dynamics requires Hill-type non- 
linearity in order to show bistability and switching; 
for the stochastic dynamics, cooperativity is not 
needed. 
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B. The repressilator 

Clearly, the dynamics of the genes does affect the sys- 
tems transients, and as such the genes can indeed have a 
profound influence on the dynamics, as we now show. For 
this we look at a gene circuit whose stationary behaviour 
is not governed by a simple fixed-point, but by a limit cy- 
cle: the repressilator. The repressilator is the three-gene 
negative- feedback loop shown in Figure 2.1; this system 
has been realized experimentally as a synthetic gene cir- 
cuit in E. coli [17] . and it has recently been the topic of 
various modeling papers, employin g b oth deterministic 
and stochastic approaches, e.g., [hH. |24|. l25l [2^| . 

The nonlinear dynamics of the repressilator in the null- 
cline space of the gates is described by the ode 




1 + vb h 



7a 



(27) 



with the equations for b and c to be obtained from the 
permutations (a — ► b,b — ► c) and (a — * c,b — > a). 
Since all parameters are assumed equal the system has 
a symmetric fixed-point, a = b = c = a. Testing the 
stability of this fixed-point the stability matrix reads as 



p/p — 



"7 




(28) 



with k = ehva h 1 /(1 + va h ) 2 . The characteristic poly- 
nomial to this matrix is given by 



(7 + A) + K 



.3 



0. 



so that the first eigenvalue is found to be 



Ai 



-(7 + k). 



The two others are given by 

A2,3 



fx fx l — 

7 + ^^. 



The condition for a Hopf-bifurcation therefore is 



2= 7 ' 



(29) 



(30) 



(31) 



(32) 



Making use of the fixed-point conditions one finds the 
relation 



(33) 



and hence the condition on the Hill-exponent h > 2 for 
the circuit in order to have a stable limit cycle. 

The stability analysis of this fixed-point can be carried 
out analytically for the full gene gate circuit, i.e. keep- 
ing both the transcription factors and the three genes as 
dynamic variables. By symmetry, in fact, the calculation 
works for a circular circuit of n genes. The calculation 




Figure 4: Top: The repressilator dynamics without gene gates 
(fixed at the nullclines of the gates) for the parameters r = 1, 
7 = 0.1, e = 0.3, r\ = 0.9, h = 3: the limit cycle is absent, 
the fixed-point is stable. Bottom: Plot of the repressilator 
dynamics for the full system with identical parameters: the 
limit cycle persists in a wider range of parameters. 



amounts to generalize eq.(26) so that 

(7 + A)"(A + X ) n + (-l) n -\sO n = ■ (34) 

with n = 3 for the repressilator. This fixed-point condi- 
tion is formally equivalent to that of the "leaky" repressi- 
lator discussed in [l4[ , for which a condition h > 4/3 was 
established. Within the full gene gate dynamics, the con- 
dition on h is thus weakened: the repressilator already 
oscillates for Hill exponent values less than two. Even 
for the case h = 3, e.g., when both the full and the re- 
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s 



Figure 5: Parameter regimes for the repressilator dynamics. 
I: stable fixed-point; II: stable fixed-point for the reduced sys- 
tem, limit cycle for the full system; III: limit cycle. Parame- 
ters are as in Figure 4. 



stricted system show oscillatory behaviour, the presence 
of the gene dynamics enlarges the oscillatory region in the 
space of protein concentrations. The stability of the limit 
cycle in the space of parameters (e, rf) is summarized in 
Figure 5. 



By contrast, the stochastic repressilator without coop- 
erativity displays a limit cycle behaviour, as shown in 
Figure 6 (top). The limit cycle appears as a symmetric 
triangle in the space of transcription factor concentra- 
tions (a, b, c). The triangle is somewhat 'fuzzy', reflecting 
the fluctuating nature of the concentrations. This fuzzi- 
ness can be reduced by increasing the space of variables 
in the system. In a recent study, the effect of an inclu- 
sion of transcription factor cooperativity (dimerization 
and higher), or an inclusion of explicit RNA transcrip- 
tion and protein translation was studied. It was found 
that all these mechanisms regularize the oscillatory be- 
haviour [25| and render the limit cycle less 'fuzzy'. Anal- 
ogous findings for circadian clocks were reported earlier 
[30l [3l| . The corresponding limit cycle for the deter- 
ministic dynamics of the reduced system is shown in the 
bottom graph. Here again a Hill coefficient h = 3 has 
been assumed. 




Figure 6: (Color online) Top: The limit cycle of the stochas- 
tic repressilator. Simulation parameters are: r = r p — 1, e 
— 0.1, i]— 10 -2 , 7 = 10 -3 . Bottom: the deterministic ver- 
sion for comparison (reduced system in region III of Figure 5, 
parameters identical to the stochastic version, with h = 3). 



IV. MULTI-INPUT GATES 

A. A rewired repressilator 

The 'stabilizing' effect due to the presence of the gene 
gates persists in the presence of multiple inputs, in fact, 
in can even be reinforced. We observed this when consid- 
ering a rewired repressilator shown in Figure 7, in which 
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with 



"1 



1 




Figure 7: (Color online) Top: the rewired repressilator: a 
positive loop is added (see arrow), so that one of the genes 
is doubly regulated. Bottom: the limit cycle of the (reduced) 
rewired repressilator circuit; the additional activation inter- 
action breaks the symmetry, as discernable in the difference 
in maximal concentrations. Simulation parameters are: r —1, 



10" 



e = 0.1, rn = 772= 10" 2 , 7 = 10-*, h = 3). 



an additional activatory loop has been added so that we 
have 



neg(c; b)\posneg(c, 6; a)\neg(a; c) 



(35) 



In the case without genes, this means that one of the 
equations, say the one for a is replaced by 



Kq 



vhb h - 1 (e + r p c h ) 
~ (1 + vb h + v p c h ) 2 



Ki 



ftc' l - 1 (rp(l + vb h ) - v p e) 



(1 + vb h + v p c h ) 2 



uhc 



h-l 



K2 = ~ 



vhb 



h-i 



K3 



(1 + vb h f 



(38) 



(39) 



(40) 



(41) 



Note that K\ can be both positive and negative. The 
characteristic polynomial reads 



(7 + A) 3 + (7 + A)kik 3 + = 



(42) 



which still has a pair of complex eigenvalues. The Hopf 
condition is given by 



87 + 27K1K3 — K0K2K3 = . 



(43) 



The analysis of the full system, genes included, is 
clearly more involved than for the repressilator due to the 
increased number of variables. We have therefore studied 
the system only numerically and compared the reduced 
and the full version, as we did for the repressilator. Our 
calcuations show that the reduced version (3 ode's for 
a,b,c) is less robust against rewiring than the gene gate 
version (7 ode's): the stability limit of the limit cycle 
regime can differ by parameter values up to one order 
of magnitude. This finding is notable since in the pres- 
ence of multiple regulations the number of gene states 
increases linearly with the number of inputs (neglecting 
still additional regulatory layers) and thus significantly 
enhances the complexity in modeling circuits with such 
elements. We close the section with Figure 7 (bottom) 
which shows the limit cyle of the rewired repressilator for 
the reduced deterministic system (h = 3). It illustrates 
that in general the presence of the additional positive 
loop breaks the (a — b—c) symmetry between concentra- 
tions. 



e + r p c" 



1 + vb h + v v c h 



7a 



This 'rewired' repressilator still has a unique fixed-point 
(a,b,c), as follows from an analysis of the fixed-point 
conditions. The stability condition can be read off, as 
before, from the stability matrix which now reads as 



p/p _ 




(36) B. A multi-input circuit related to developmental 

regulation 

In this final subsection we address a second example 
of a multi-input gate. It consists of a bistable switch 
built from two repressing gates which is placed under ad- 
ditional control by an activating input. Such motifs oc- 
cur both in transcriptional regulation [l|| , but they have 

(37) also been proposed recently to play a role in morphogen 
concentration-dependent cellular development [l8j; our 
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example is motivated by the latter case. The circuit dy- 
namics is governed by the following ode's (neglecting the 
gene gate dynamics since we are concerned with fixed- 
point dynamics only) 



b = 



ra 



- vc m + i> ac a r ' 
e c + ra n 



1 



b< 



Vaba'' 



- 76 
7c 



(44) 
(45) 



where m, n, I are the different Hill exponents. If the acti- 
vating variable a = 0, the system is the standard bistable 
switch, albeit asymmetric with respect to the parameters 
and nonlinearities, and it is this asymmetry which plays 
an important role - in ref. [18j, the supposed Hill coeffi- 
cients have values of 3 and 6, respectively. 

The effect of the variable a has on the dynamics is eas- 
ily understood. To simplify matters, we neglect a in the 
first equation and look at an asymmetric wiring. It actu- 
ally does not matter whether we allow a to control one or 
both transcription factors b, c as long as a interacts with 
both in the same way and not via a different nonlinear- 
ity: the main symmetry-breaking effect is contained in 
the difference between the Hill coefficients controlling b 
and c. 

Supposing further that we increase the concentration 
of a to levels where it dominates the concentration b so 
that we have for the fixed-point in c 



1 e c + ra n 
Co = - 



Vab 

7 1 + v a ba n 7 



(46) 



Thus, the fixed-point concentration of the repressing vari- 
able Co is locked to that of a and approaches an asymp- 
totically constant value. Correspondingly, this brings the 
fixed-point level of b down and under firm control of a: 
the system ceases to be bistable, and locks into a stable 
state under control of a. The possible relevance of this 
mechanism for a transcriptional circuit in development is 
evident: increasing a can force the system to switch in a 
concentration-dependent way. 

In the nonlinear dynamics case, this switch is there- 
fore brought about by the vanishing of a fixed-point; 
again, this situation is different in the stochastic setting. 
For comparison, Figures 8 and 9 show our results of the 
stochastic simulations for the circuit 
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Figure 8: (Color online) Top: the bistable switch under ex- 
ternal control by a. Bottom: the system starts at zero con- 
centrations of both proteins and enters the state with higher 
stability, as given by higher production rates. The signal c 
dominates widely over a, b which are indistinguishable from 
the baseline. 



parameter range, in which bursts in concentration c can 
occur at random times within a wide time interval (see 
the concentration peak at around 55.000 a.u.), and are 
finally controlled by a. In Figure 9 (bottom) the system 
has switched to a dominant concentration of b and the 
concentration of the previously dominant transcription 
factor c is now fully controlled by a. Note the difference 
in concentration levels of all proteins in the figures. 



null (a) \posneg{a, b; c) \posneg(a, c; b) (47) 

without any cooperative nonlinearity, as for the repres- 
silator. The progression of dynamic behaviours in Fig- 
ure 8 bottom to Figure 9 top and bottom is controlled 
by the average concentration level of a, which increases 
from one figure to the next by one order of magnitude 
since the transcription rate is increased by this factor. 
In Figure 8 (bottom) the switch enters the more stable 
of the two states; in Figure 9 (top) the additional input 
a makes the concentration levels b and c compete with 
each other. This behaviour is observed within a large 



V. DISCUSSION AND OUTLOOK 

In conclusion we have proposed a minimal model de- 
scription for gene regulatory networks based on the no- 
tion of the gene gate, first proposed in ref. [24]. We 
studied the dynamics of simple gene networks in both a 
mean-field and a stochastic version, with characteristi- 
cally different results: 

• If the system dynamics is stable fixed-point only, a 
reduced deterministic description ignoring the de- 
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scription fluctuation effects induced by the genes 
(promoters) might affect fixed-point locations 32 1, 
or the stability of bistable switches [H, [H, HH, [36 1 . 
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• If the system displays a limit cyle, the gene gates 
are relevant, as is any other additional regulatory 
layer to determine the parameter range of oscilla- 
tions. In general the limit cycle regimes depends on 
the whole set of parameter values, Hill coefficient 
included. In particular this means that in multi- 
input regulations in which additional gene states 
have to be accounted for, the parameter space can 
extend significantly. 

• If the system dynamics is fixed-point, the stochas- 
tic version obeys this without any need for cooper- 
ative effects. The same holds true for limit cyle be- 
haviour. Additional regulatory layers also enlarge 
the phase space but in a trivial way. By contrast, 
they affect oscillatory behaviour by regularizing the 
oscillations. 

In our view these results have interesting consequences 
on the philosphy of modeling gene regulatory networks in 
suggesting a different coarse-approach. Computational 
models of large networks can be built by abstracting 
away all regulatory layers to a level where the remaining 
network can still faithfully represent the system char- 
acteristics. Network motifs that have a more sensitive 
dynamic behaviour - like limit cycles, as shown here - 
are more sensitive to modeling assumptions. Finally, we 
remark that in view of our results, modeling attempts 
combining deterministic and stochastic aspects should 
be considered with care [371. 
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Figure 9: (Color online) Top: a change of the activator tran- 
scription rate by one order of magnitude makes both proteins 
compete; note the concentration overshoots of the previously 
stable protein. Bottom: A further increase of the activator 
transcription rate makes the system switch between the two 
states. Simulation parameters are (bottom): r =1, e a = 10 -3 , 
£ 6 = 0.1, e c = l(r s , r, ab = 2-HT 3 , r, ac = 2-HT 2 , t? = 2-HT 1 , 
7a = 5-HT 3 , 7 = 2 -lO" 4 . 



grcc of freedom of the gates is sufficient in the sense 
that the fixed-point is not altered by the presence 
of the genes. But if this is the case, the latter are 
indeed 'irrelevant'. In order to represent faithfully 
the fixed-point structure of the network, a Hill-type 
nonlinearity may be needed (like for the bistable 
switch circuit). However, within a stochastic de- 
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VI. APPENDIX: SIMULATIONS IN 
STOCHASTIC 7T-CALCULUS 

The Gillespie algorithm can, of course, be implemented 
in various different programming languages. What then 
are the main ideas and advantages of the 7r-calculus? 

The 7r-calculus is a formal system in which each compu- 
tation is represented by a communication over input and 
output channels. The communicating objects are called 
'processes'. Computation by communication within pi- 
calculus can be understood as an alternative to, e.g., 
functional computation as realized in the A-calculus. The 
7r-calculus is Turing complete: it can therefore realize any 
possible computation [22j. 
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For our application, the calculus allows to represent 
each gene gate by a computational process 

gate{x; y) (48) 

with its corresponding input(s) x and output(s) y; e.g. 
the repressing gate of Figure 1.2 is written as neg{a;b) 
where the input channel a represents the repression of 
transcription by transcription factor a, and b is the cor- 
responding output. All other reactions, like e.g. the 
degradation process of b, are bound to this process and 
contained in its definition. 

The scheduling of inputs and outputs on a gate are 
calculated in the usual fashion by the standard Gillespie 
algorithm, as adapted to the 7r-calculus [24l l25l l26l|. 

One main technical advantage of the calculus is, in fact, 
that its syntax and semantics are perfectly adapted to a 
'compositional' build-up of the transcriptional networks. 



In our context this permits to express (and compute!) 
a composed circuit, like the repressilator, by a parallel 
process 

neg(c; b)\neg(b; a) \neg(a; c) (49) 

The second main advantage (although not exploited for 
the small systems studied here) is that it can reduce the 
computational complexity of a system of n kinetic reac- 
tions, which is of order n 2 , to linear order. The interested 
reader is referred to refs. [UHII for more details, writ- 
ten in a way accessible to a physics-trained audience. The 
simulation results presented here were obtained with the 
public domain software SPIM, downloadable with doc- 
umentation and examples (27[. The details of the im- 
plementation of the Gillespie algorithm in the dedicated 
software SPIM are discussed in the Supplementary Ma- 

terials of HH. 



[19] 



A. Goldbeter, Biochemical Oscillations and Cellu- [20 
lar Rhythms. Cambridge University Press, Cambridge, [21 
United Kingdom (1996) 

CP. Fall, E.S. Marland, J.M. Wagner and J.J. Tyson [22 
(eds.) Computational Cell Biology. Springer, Heidelberg, 
Germany (2002) [23 
N. van Kampen, Stochastic processes in physics and 
chemistry, North-Holland, Amsterdam, The Netherlands [24 
(2007). 

D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977) [25 
J. M. G. Vilar, H. Y. Kueh, N. Baarkai and S. Leibler, 
Proc. Natl. Acad. Sci. USA 99, 5988 (2002) [26 
P. S. Swain, M. B. Elowitz and E. D. Siggia, Proc. Natl. 
Acad. Sci. USA 99, 12795 (2002) [27 
M. B. Elowitz, A. J. Levine, E. D. Siggia and P. S. Swain, [28 
Science 297, 1183 (2002) 
J. Paulsson, Nature 427, 415 (2004) [29 
J. Paulsson, Physics of Life Rev. 2, 157 (2005) 
H. de Jong, J. Geiselmann, C. Hernandez and M. Page, [30 
Bioinformatics 19, 336 (2003) 

L. Glass and S. A. Kauffman, J. Theor. Biol. 39, 103 [31 
(1973) 

R. Thomas, J. Theor. Biol. 42, 565 (1973) [32 
R. Bundschuh, F. Hayot and C. Jayaprakash, Biophys. 
J. 84, 1606 (2003) [33 
S. Miiller, J. Hofbauer, L. Endler, C. Flamm, S. Widder 
and P. Schuster, J. Math. Biol. 53, 905 (2006) [34 
S. Widder, J. Schicho and P. Schuster, J. Theor. Biol. 
246, 395 (2007) [35 
U. Alon An Introduction to Systems Biology, CRC Chap- 
man & Hall, London, United Kingdom (2006) [36 
M.B. Elowitz and S. Leibler, Nature 403, 335 (2000) 
Y. Saka and J. C. Smith, BMC Developmental Biology [37 
7, 47 (2007) 

M. Madan Babu and S.A. Teichmann, Nucl. Acids Res. 
31, 1234 (2003) 



J. N. Weiss, The FASEB Journ. 11, 835 (1997) 

J.L. Cherry and F.R. Adler, J. Fheor. Biol. 203, 117 

(2000) 

R. Milner, Communicating and mobile systems: the ty- 
calculus, Cambridge University Press (1999) 

C. Priami, A. Regev, W. Silverman and E. Shapiro, Inf. 
Proc. Lett. 80, 25 (2001) 

R. Blossey, L. Cardelli and A. Phillips, T. Comp. Sys. 
Biology IV, 99 (2006) 

R. Blossey, L. Cardelli and A. Phillips, HFSP Journal 2, 
17 (2008) 

A. Phillips and L. Cardelli, Comp. Meth. Sys. Biol. 4695, 

184 (2007) 

http: / /research. microso ft.com/^aphillip/l 

A. Lipshtat, A. Loinger, N. Q. Balaban and O. Biham 

Phys. Rev. Lett. 96, 188101 (2006) 

A. Loinger and O. Biham, Phys. Rev. E 76, 051917 
(2007) 

D. Gonze, J. Halloy and A. Goldbeter, Proc. Natl. Acad. 
Sci. USA 99, 673 (2002) 

D. Gonze, J. Halloy and P. Gaspard, J. Chem. Phys. 116, 
10997 (2002) 

J. Paulsson, O. G. Berg and M. Ehrenberg, Proc. Natl. 
Acad. Sci. 97, 7148 (2004) 

P. B. Warren and P. R. ten Wolde, Phys. Rev. Lett. 92, 
128101 (2004) 

R. J. Allen, P. B. Warren and P. R. ten Wolde, Phys. 
Rev. Lett. 94, 1804 (2005) 

A. M. Walczak, J. Onuchic and P. G. Wolynes, Proc. 
Natl. Acad. Sci. 102, 18926 (2005) 

M. J. Morelli, S. Tanase-Nicola, R. J. Allen and P. R. ten 
Wolde, Biophysical J. 94, 3413 (2008) 
M. Scott, T. Hwa and B. Ingalls, Proc. Natl. Acad. Sci. 
USA 104, 7402 (2007) 



